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Abstract 

In this paper, we describe a new class of fast solvers for separable elliptic partial differential 
equations in cylindrical coordinates {r,9,z) with free-space radiation conditions. By combining 
integral equation methods in the radial variable r with Fourier methods in and z, we show 
that high-order accuracy can be achieved in both the governing potential and its derivatives. 
A weak singularity arises in the Fourier transform with respect to z that is handled with special 
purpose quadratures. We show how these solvers can be applied to the evaluation of the Coulomb 
collision operator in kinetic models of ionized gases. 
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1. Introduction 

A variety of problems in computational physics require the solution of the Poisson and bihar- 
monic equations in cylindrical coordinate systems, particularly when the source distribution (the 
right-hand side) is axisymmetric or involves only a few azimuthal modes. The present paper was 
motivated by the need to compute the Coulomb collision operator C{fa,fb) in kinetic simulations 
of the Boltzmann-Fokker-Planck equation fl H [H [Hill] : 



d.fa +V-Vfa + —{E + VXB)-dyfa^y C(Ufb). (1) 

ma ^ 

Here, fa{x, v, f) denotes the state of an ionized gas for plasma species a and the index b runs over 
all species present. In the Fokker-Planck-Landau formalism 12611 . 



C{fa,ft)-yatd.- r §(V-V')(^^/.(V')-/»^^^)^V' (2) 

J \ ma nib I 



where 

1 (vi - v'i){vj - v'p 
|v-v'| |v-v'. 
An alternative representation makes use of the Rosenbluth potentials Js^]: 



- = '^'VnrTT^ ,,' ..,,3 ' ■ (3) 
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C{fa,fb)=—dy 



■ (fad,d,Gb) -2(1 + ^) fMh 

ml 



where 



and 



Hb(v) 



Gb(v) 



/i 
/ 



-Jb(v')dv' 
It' 



or 



AHh 



v-v'\Mv')dv' or A-Gb 



-4nfb 



-Snfb 



(4) 



(5) 



(6) 



Note that four derivatives of Gh are required in (|4}, while Gb itself satisfies the inhomogeneous 
biharmonic equation (|6]l. Thus, direct discretization of the partial differential equation, followed 
by evaluation of the collision operator via (|4) would require eight steps of numerical differentia- 
tion, with significant loss of accuracy. 

It is natural, therefore, to consider alternative methods with the dual goals of achieving high 
order accuracy and minimizing the condition number of the solution process. Because of the 
design of magnetic confinement devices for plasmas, it is also important to be able to construct 
numerical methods in cylindrical coordinate systems, since the distribution functions fbiv) are 
often axisymmetric or involve only a few azimuthal modes. 

There is, of course, a substantial literature on computing Coulomb collisions and on solving 
elli ptic pa rtial differential equations in cylindrical coordinates. We refer the reader to llsl ItLIiiI 



Mm 



M^MMM,^ for some methods in current use in plasma physics. For a discussion 

3C 



of relativistic effects, see la]. Most closely related to our approach are the methods of 1115 



and 12^ 23, 24, 32]. The first two are fast and achieve high order ("spectral") accuracy, but use 
Fourier methods in Cartesian coordinates and do not address the axisymmetric (or low azimuthal 
mode) case. The latter rely on separation of variables in spherical coordinates, for which the 
axisymmetric case leads naturally to a representation involving Legendre polynomials and the 
general case to a representation involving associated Legendre functions. 

In the numerical analysis literature, most solvers based on cylindrical coordinates tend to con- 
cern themselves with periodic (in z) or finite domain boundary conditions rather than free-space 
boundary conditions (see, for example ^ |25|). Here, we develop a method for computing the 
Rosenbluth potentials using separation of variables and a mix of integral equation and Fourier 
analysis techniques. We show that free-space (radiation) conditions can be imposed in a straight- 
forward manner and that high order accuracy can be achieved in all derivatives with minimal loss 
of precision. The solver requires 0{N log A^) work, where is the number of grid points used to 
sample the distribution function. 

Finally, we should make a remark about notation. The collision operator and the Rosenbluth po- 
tentials in (|5]l,(|6]l are defined in velocity variables, for which we will use the standard cylindrical 
coordinates (r, 6, z) for v. In the context of plasma physics, r = | Vj^|, where Iv^l is the magnitude 
of the component of the velocity perpendicular to the magnetic field, 6 is the gyrophase angle, 
and z - Vft is the component of the velocity field parallel to the magnetic field. The problem is 
purely axisymmetric when the velocity field is independent of the gyrophase angle. 
One disadvantage of our solver is that we can be adaptive in the r direction, but not in the z or 
directions, since we use spectral discretizations in the latter variables. For fully adaptive three- 
dimensional calculations, one could employ fast multipole-accelerated integral equation solvers, 
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as described in lfl4l EtIi . These methods directly compute the convolution of the data fh(v) 
with the free-space Green's function. In the axisymmetric case, one could use an axisymmetric 



version of the fast multipole method 113411 . The constant, however, is larger for these schemes 
than for methods based on separation of variables, and we limit our attention to methods that 
rely on a tensor product mesh in r, 6 and z, which is adequate for most current simulations of the 
Boltzmann-Fokker-Planck equation ([TJ. 

2. The Poisson equation in cylindrical coordinates 

In order to compute the Rosenbluth potential we must solve the Poisson equation in free 
space 

Au{v) = /(V). 
In cylindrical coordinates v - (r, 0, z), we have 

Urrir, e, z) + -Urir, 8, z) + \ugg{r, 6, z) + u-Ar, 0, z) = f(r, 0, z), (7) 
and we assume that / is identically zero outside the region 

n^{{r,9,z):0<r<R, -A < z < A, 0<6<2n}. 
Since u and / are periodic in 0, we represent them as Fourier series: 

oo 

u(r,0,z) = 2 u^"\r,z)e'"'' (8) 

— CO 
CO 

/(r,0,z)= /"V-zy"" (9) 

«=-oo 

The derivatives in this representation will be written as 

CO CO 

n=— CO /7— — CO 

OO CO 

u,^{r, 0, z) = ^y"' "eeir, e,z)^ ^ (-nW^Hr, z)e'"» 

n=-oo n=-oo 

Substituting into d?) and equating terms corresponding to the nth azimuthal mode, we obtain: 

M<';)(r,z) + -M("'(r,z) - + u'-"\r,z) = f"\r,z). 

r 

For each mode, we now have a partial differential equation (PDE) in the two variables r and z 
which we need to solve on the rectangular domain 

^r, = \(r,z) ■.0<r<R, -A<z<A}. 
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Let us now take the Fourier transform of the equation in the z direction. That is we write 



1 

„(")(r,z) = — e'''~u^"\r,K)dK (10) 
e-''^-u^"\r,z)dz 

OO 

1 r°° 

f"\r,z)^- e'''f"\r,K)dK (11) 

/-*oo 

f"\r,K)^ e-''"f"\r,z)dz 



In the Fourier transform domain, the PDE becomes an ordinary differential equation (ODE), 
where k (as well as n) is now fixed: 

u^l\r, k) + ^u'-;\r, k)-^-^+k^^ u^"\r, k) = f"\r, k). (12) 

To simplify notation (when the context is clear), we will write u(r) instead of u^"\r, k) and u'{r) 
instead of M|"'(r, k) to denote the derivative when discussing the solution of the ODE. 
The equation ( fT2b is an inhomogeneous modified Bessel equation 12]. In the homogeneous case, 
the equation has two linearly independent solutions, namely IniMr) and Kn{\K\r), the modified 
Bessel functions of order n. The function I„{\K\r) is regular at the origin, and grows exponentially 
as r — > oo, while K„{\K\r) is logarithmically singular at the origin, but decays exponentially fast 
as r — > oo. 



2.1. Boundary conditions for the modified Bessel equation 

In order to have a properly posed ODE, we seek two boundary conditions, one at r = and one 
?Ar - R, beyond which the equation is homogeneous. For the n - Q mode, the condition 

m'(0) = uf\Q,K) = 

ensures regularity at the origin, while for modes n + Q 

m(0) = u^"\Q, k)^Q 

is necessary. This is easily seen from taking the limit of the equation (fT2l i as r — » 0. 

Since we are seeking to solve the Poisson equation in free space, our ODE is actually posed on 

the half line [0, oo], with the radiation condition that the solution decay at infinity. This can be 

accounted for exactly in terms of a suitable boundary condition r - R. To see this, note that 

for r > R the solution must be proportional to Kn{\K\r), since In{K\r) grows without bound. That 

is, 

u{r)^C„,,-K„{\K\r) fov r > R, 

where C„ ,f is an unknown constant. The solution on [0, R] and its derivative must match this 
solution at r = /?, so that 

u'{R) = Cn,M ■ K(\k\R). 
4 



Eliminating the constant C„,,f we obtain the exact "radiation" boundary condition: 

In summary, the ODE boundary value problem we must solve is (fTSl i. subject to the boundary 
conditions: 

m'(0) = n = 

m(0) = « 5t (14) 

~ rr^^TTn^" = o. (i5) 
kl ■ ■'^//kl^) 

In broad terms, this completes the description of the Poisson solver, which proceeds in four steps. 



Informal description of fast Poisson solver 

1. Expand the right hand size /(r, 0, z) as a Fourier series in 0, in order to get f^"\r, z), 

2. Compute the Fourier transform of f^"\r, z) in the z direction to get f^"\r, k), 

3. Solve the ODE ( fT2l i for each k and n to obtain M^"'(r, k), 

4. Compute the inverse Fourier transform of M*"*(r, k) to get M*"*(r, z), 

5. Sum the Fourier series in to get the final solution u(r, 0, z). 



We will rely on fairly standard methods for all of the above, except Steps 3 and 4. For Step 3, 
we use an analytic solution based on knowledge of the underlying Green's function for the ODE 
and accelerated by a simple "sweeping" algorithm. Step 4 will require some care, since it is 
straightforward to show that uS"\r,K) is logarithmically singular as /c — > for n = and has a 
singularity of the order a:^'"' log k for n i^Q. 



3. Discretization and solution 



We assume /(r, 0, z) is given on a tensor product grid with Ne equispaced points in the direction 
on [0, 2;7r], A^^ equispaced points in the z direction on [-A,A], and A^,. points in the r direction 
on [0,/J]. We divide [0,/?] into Ni intervals with interval endpoints - 0,RuR2, ■ ■ ■ ,Rn, - R- 
We use a f th order (scaled) Chebyshev grid on each, so that A^, = Ni P. We will denote by 
{rj \ j - 1, ;A^,.) the grid points in increasing order When the particular interval m (1 < m < Nj) is 
of interest, the pth grid point on that interval (1 < p < P) is rj = r^m-i) p+p- 
The discretized data will be denoted by 

Mn, 0n,Zk) = f(rj, 0„,Zk) for < 7 < A?,, < « < A^e, < ^ < N,, 
3.1. Step 1: Transformation in 

We use the fast Fourier transform (FFT) to compute fl"\rj, Zk), the discretized version of f^"\r, z). 

f}"\rj, Zk)^^Yj ^"^"'//.(O- ^k) « f"\rj, Zk). (16) 

" 1=0 
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It should be noted that, if /(r, 6, z) is n-times differentiable, then the series (|9) trancated after 
terms has an error of the order 

^(a^) ^^^^ 

If / is infinitely differentiable, then the error goes to zero faster than any finite power of 1 jN. 
Schemes with this property are often referred to as having spectral accuracy. Moreover, the 



trapezoidal rule approximations of the series coefficients in ( fT6b converge at the same rate 116 



3511. 



3.2. Step 2: Transformation in z 

Since f{r, 9, z) and f"\r, z) are compactly supported, we need to compute the finite integral 

f"\r,K)^ r e-''' f"\r,z)dz. 



Letting h, - ^ and zi - Ih^, the trapezoidal rule yields: 



2A 

l=-NJ2 



Z^) = f Z ^-^'^^rV.z,) ^ />' (n, ) . (18) 



This is computable using the FFT, and yields the values of the Fourier transform at equally spaced 
points of step size n/A in the k domain. A few remarks are in order: 

• The ratio ^ determines the range of frequencies that are resolved. If ^ increases 
decreases), higher frequency modes of the data are computed. 

• We will assume that, to precision e, f{r, k) is supported on the interval 1-/Cmax, fmax]. For a 
given A, A^; must be chosen sufficiently large that nN^JilA) > k,^^^. (This is simply asking 
that the grid in z be fine enough to resolve the data.) 

• Increasing A^- and A simultaneously so that N~/A remains fixed leaves the range of k un- 
changed, but increases the number of sample points where is computed in the range 
t-7rA^;/(2A),7rA^./(2A)]. 

• The trapezoidal approximation (ITST l is spectrally accurate, since the integrand and all its 
derivatives are assumed to have vanished by the time z - ±A. 

3.3. Step 3: Solving the modified Bessel equation 

We turn now to the solution of the modified Bessel equation (fT2l l. subject to the boundary condi- 
tions (ITST l for K ^ 0. (As noted above, the equation has a weakly singular solution at a- = 0. Our 
quadrature rule for computing the inverse Fourier transform in section l34l will avoid the origin 
when integrating along the k axis.) 

One possible approach to solving the equation is to use a spectral integration-based ODE solver 
1 1711 that represents the second derivative as a Chebyshev series: 



M"(r) = ^a,r,(r). 



k=0 
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Multiplying the equation (fT2l) by r~ and systematic use of the following two identities for Cheby- 
shev polynomials 



Tn{r)dr + C = ^ T„^i(r) - — ^— r„-i(r) 
2(« + 1) 2{n - 1) 

r„+,(r) + r„_i(r) 
rT„(r) = — 



2 

yields a banded linear system (of bandwidth 7) to which are appended two dense rows that 
correspond to the imposition of the desired boundary conditions. Such a system can be solved 
in linear time by careful Gaussian elimination, achieving spectral accuracy. For non-singular 
ODEs, this linear system can be viewed as the discretization of a second-kind integral equation 
for the unknown second derivative, and thus as a well-conditioned formulation of the problem. 
Unfortunately, in our case, the differential operator is singular at the origin. As a result, the 
integral equation is not of the second kind and the approach becomes ill-conditioned for fine 
grids, with the attendant loss of precision. 

An alternative strategy is to use the fact that our ODE is classical and well studied, with a known 
Green's function G"{r, s). We can, therefore, write down the exact solution as a convolution: 



Jo 



u'"\r,K)^ I G:{r,s)f{s)ds (19) 
Jo 



where 



^ I„(Kr)K„iKs)/W(s) if r<s , 
GJr, i) = < where 
[ K„(Kr)I„(Ks)/W(s) if s<r 

Wis) = k(I;,(ks)K„(ks) - K'„(ks)Uks)) = -- 

s 

This choice of Green's function correctly imposes the regularity condition at the origin and the 
radiation condition at infinity. In this formulation, there is no need to solve a linear system - one 
needs only to evaluate the integral in (1% . Naive implementation of this formula would require 
0(Nf) work. Because of the structure of the Green's function, however, there is a simple 0(Nr) 
solver based on the observation that 



u^"\r,K)^K„(Kr) I„(Ks)f(s)/W(s)ds + 
Jo 

UKr) f K„(Ks)f(s)/W(s)ds . (20) 



The only source of error comes from the quadrature approximation of the preceding integrals. 
Derivatives of the solution are also obtained analytically. For example. 



u'-;\r,K)^KK'„(Kr) f UKs)f(s)/W(s)ds + 
Jo 

Kl',{Kr)J^ Kn{Ks)f(s)IW(s)ds. (21) 



There are some implementation issues in using (120) . having to do with scaling and quadrature 
due to the fast growth/decay of Bessel functions for increasing n and r. In particular, when r lies 
in the mth interval denoted by [7?„,_i,7?„,], we write 

u^"\r,K)^ ^'f'''\ f I„{Ks)K„{KR„,)f{s)/W{s)ds + 
k„(kK,„) Jo 

UKr)K„iKRJ "\\ fis)/Wis)ds. (22) 

J,. K„(kR„,) 

3.4. Step 4: Computing the inverse Fourier transform 

We now need to compute the inverse Fourier transform of u^"\r, k) to recover u^"Hf, z), according 
to ( fTOb . Since u^"\r,K) is compactly supported to the desired precision on [-A:niax, 'Cmax], we 
actually need to compute 



1 rs^= 

„W(r,7)^_ e""i/"\r,K)dK (23) 

27T 



where (as discussed in section [T2l i nN^I{2A) > K„,„y. A complication is that ii^"\r,K) has a 
logarithmic singularity at k - 0. 

Fortunately, in the last decade or so, a variety of quadrature rules have been developed that rely 
on slight modifications of the trapezoidal rule, yield high-order accuracy, and still permit the use 
of the FFT. Two such schemes are the end-point corrected trapezoidal rule due to Kapur and 
Rokhlin fl2lll and the hybrid Gauss-trapezoidal rule due to Alpert 130. We will make use of the 
latter 

Theorem 1. (modified from Let f(K) be a compactly supported function on [—Kmax,Kniax\ 
which is smooth away from the origin and takes the form 

f{K) = Sl{K)log{\K\)+ S2{k) 

in a neighborhood of the origin, where si and si are smooth functions. Let 

/(/)= ^ f{K)dK 

and let h — Then, for every integer m > and every > 2m, there exist weights wij„ and 

nodes K^m such that 

NJ2 in 

hif)^h £ fik^+Y.y^i.nfM (24) 

k=-NJ2 l=-m 
\k\>m 

satisfies 

ih(f) ^ Kf) + 0(h'") . 

In other words, the hybrid Gauss-trapezoidal rule achieves m''' order accuracy by replacing the 
2m trapezoidal nodes nearest the origin with specially located nodes (and weights). The paper 
istl provides tables of these nodes for orders 2-16 (and the corresponding ones for a variety of 
other singularities as well). 
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In the present context, therefore, we will compute the integral (l23T l using the formula (l24l : 



, NJ2 m 

u"\r,z) Tj ^""«'"*('-'^*) + Yj ^i,'ne'""''u^"\r,K„n) (25) 

k=-NJ2 l=-m 
\k\>m 



u^"\r, z) on our grid Zj - ^ j, we have: 



with mesh spacing h - 2nN-l{2A) ^ N~, - jk, and wi,„,ki^,„ taken from yD. Evaluating 

2A 

2A 



, NJ2 m 



k=-N,/2 l=-n 
\k\>m 

The first term is straightforward to compute with the (inverse) FFT, requiring N^logN^ opera- 
tions. The second term can be computed directly using 0(Njn) operations, where m is the order 
of the quadrature rule. For m sufficiently large, the sums can be computed simultaneously using 



the non-uniform FFT (see \ and the more recent review 1 18 1 ) 



Remark 1. The quadrature rule ( 1251 ) determines the discrete values of the continuous Fourier 
transform variable k where u^"\r, k) needs to be sampled. The number of such points is 0(N,+m). 
This, in turn, tells us where f^"\r, k) is needed. The values at the regular nodes are obtained with 
the FFT, as discussed in section 13.21 The values at the irregular nodes Kim can be computed 
directly or using the non-uniform FFT. 

Remark 2. (Oversampling) . In practice, there is one more issue which needs to be addressed. 
In the integral f l2il ), z is bounded by A, so that the most oscillatory integrand is e""^u^"\r, k). It 
is easy to see that there are a maximum of N-12 periods of the function e"^'* over the interval 
of integration {-nN,l(2A),jiNJ{2A')\. The maximum for the function uf-"\r,K) is similar Thus, 
the trapezoidal rule with points yields only one point per wavelength for the most oscillatory 
argument, in violation of the Shannon sampling theorem. We, therefore, oversample the integrand 
by a factor ofrj, by setting N', — rjN^. (As discussed in section [O] we must simultaneously set 
A' — TjA in computing the forward transform.) 

Setting T] — I yields exponentially small errors near z—0, but (9(1) errors at z — A. Setting 
Tj >2 ensures convergence for z in the entire range [—A,A\, with exponential improvement as rj 
increases. Setting vj — A is sufficient for double precision accuracy for N, > 16, assuming the 
function is bandlimited to machine precision at K„,ax — nN-l(2A). 

3.5. Step 5: Sum the Fourier in 6 to obtain the full solution 

This is completely straightforward. As in Step 1, we may use the fast Fourier transform (FTT) to 
compute u{rj,6,Zk) at equispaced points 6i — 

m(0, 6u Zk)^ — J] e^t"'u'^"\rj, Zk) ■ (26) 

;7=0 
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3.6. Computing derivatives of the solution 

One useful feature of spectral solvers is that derivative are straightforward to compute with high 
order accuracy. 

1. First and second r-derivatives: Our ODE solver returns both the solution u^"\r,K) and 

(n) (n) 

its derivatives u). {r,K) and Urr (r,K) on our grid. Thus, we can compute Ur(r,6,z) and 
Urrir, 9, z) by the same technique as for u{r, 9, z): evaluating the inverse z-Fourier transform 
and the 9 Fourier series for m*"' and mJ."', respectively. 

2. z-derivatives: In the present paper, z derivatives are obtained through multipUcation by iK 
in the inverse Fourier transform step: 

1 Cim 1 /^oo 

„(«)(r,z) = — e'''u^"\r,K)dK ^ = — {iK)"'e'''u^"\r,K)dK 

and the quadrature rule described above for logarithmic singularities. 

3. 9-derivatives: In the present paper, we also compute 0-derivatives spectrally, by differenti- 
ating the Fourier series: 

u{r, 9,z)^Yj ^ ^' ^) = Z iinru^"\r, z)e'"' , 

n=—oo n=—oo 

using the FFT. 

4. Mixed derivatives: Since derivatives in r, 9, z are computed at independent steps of the 
algorithm, they are easily combined. For example, if we want to compute Urzz, we start 
with ii^^\r, k) and compute the inverse z Fourier transform on the function -K^u^"\r, k) to 
get Urzzir, z), followed by evaluating the Fourier series in the 9 direction via the FFT. 

Remark 3. One can easily obtain z derivatives without numerical differentiation, once u, Ur, u,-,- 
and Ugg are known. The original PDE becomes a second order ODE in z, and the method of 
spectral integration U^l can be applied directly. We have not implemented this option, since the 
condition number of Fourier differentiation is only 0{N), so that with 1000 points in z (or 1000 
azimuthal modes), one can still obtain at least 10 digits of accuracy in double precision. 

4. The biharmonic equation in cylindrical coordinates 

For the Rosenbluth potential Gb,vje must solve the biharmonic equation in free space 

A\{v) = /(V). 

In cyUndrical coordinates v = (r, 9, z), after Fourier transformation in z and 9, we obtain the 
fourth order Bessel-type ODE: 



) 

n'^ - 4«2 2kV 



+ /m<'"=/"> (27) 
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4.1. Boundary conditions for the fourth order Bessel-type equation 

Since (l27l i is fourth order, we need four boundary conditions to have a properly posed ODE. We 
impose two at r = and two at r = /?, beyond which the equation is homogeneous. To ensure 
regularity at the origin, for the n = mode, it is sufficient to impose 

uf\Q,K) = ufJ,iO,K) = 0. 

For the n = 1 mode, we set 

a'-^\0,K) = u';]}{0,k) = 0, 

and for modes n > 2, we set 

M'">(o,/^) = M;.'"(o,/f) = o. 

These conditions are easily derived by taking the limit of the equation ( fTSl i as r — » and the fact 
that the null space of the differential operator is spanned by 

{I„(\K\r),ri;,(\K\r),K,Mr),rK'„(\K\r)}. 

It remains to determine a radiation condition at r = R, so that the derivative of the solution is 
bounded at infinity. As for the Poisson equation, we proceed by observing that the solution u^"\r) 
for r > R must take the form 

u(r) = Cl,K„(\K\r) + ClM\K\r) , 

where C,'^, C,^^ are unknown constants. This follows since the derivatives of IiMr) and rI',{K\r) 
grow without bound. The solution on [0, R] and its derivative must match this solution at r = R, 
so that 

Cl, ■ K„(\k\R) + Cl, ■ rKUMR), 
C\, ■ kK',MV^) + Cl^ ■ (K(\k\R) + xrK'ilKlR)), 

■ k^K^\k\R) + Cl^ ■ {1kK\:{\k\R) + K^rK'^'WKm, 

■ i?K'{\k\R) + Cl, ■ (3k^K'„"(\k\R) + k\K':"{\k\R)). 

EUminating the constants, we obtain two exact "radiation" boundary conditions to be imposed 
on the combination of u^"\R) and its first three derivatives. The formula is complex and omitted, 
since we won't use it. We will instead use an exact solution based on the Green's function. 

4.2. Discretization and solution 

The solution of the biharmonic equation is analogous to that of the Poisson equation, so we just 
highlight the differences. 

After separation of variables, we need to solve a fourth order Bessel type equation. As before, 
we could proceed by expanding the highest derivative in a Chebyshev series and integrating, but 
the resulting linear system again loses precision because of the singular nature of the differential 
operator at the origin. (The loss is, in fact, much more severe than for the second order (Poisson) 
equation.) 

Alternatively, we can construct the Green's function for the ODE using the linearly independent 



u'"\R) = 

u^;\R) - 
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fundamental solutions Inixr), rI[j{Kr), K„(Kr), rKI,(Kr), imposing the regularity condition at r = 
and the decay condition as r — > oo; 

[I„{Kr)sK;,(Ks) + rr„(Kr)K„(Ks)] /Wis), r<s 2k 

with W(s) = 

[K„(Kr)sr„{Ks) + rK'„(Kr)I„(Ks)] /Wis), r>s s 

The solution involves computing four integrals (instead of two). The sweeping method is virtu- 
ally the same as that used for the Poisson equation. 

u^"\r) = K„iKr) [ sI'„iKs)fis)/Wis)ds + rK'„iKr) f I„iKs)fis)/Wis)ds 
Jo Jo 

sK'„iKs)fis)/Wis)ds + rI'„iKr) J K„iKs)fis)/Wis)ds . 

The Fourier transform of the solution in z has a more severe singularity in the biharmonic case, 
due to the fact that the free-space Green's function does not decay. Fortunately, however, we are 
only interested in second derivatives of the biharmonic potential, and they have only logarithmic 
singularities, so our special-purpose quadratures from section 13.41 yield the desired accuracy. 
(More elaborate methods involving singularity subtraction could be developed if one wanted the 
biharmonic potential or its first derivatives.) 




5. The collision operator 



Now that we've described how to solve the Poisson and biharmonic equations Q and we 
turn our attention to the collision operator (|4). If we express all the derivative terms in cylindrical 
coordinates, the axisymmetric collision operator becomes: 



cif,f) = ^ \c„if,f) - 2(1 + ^)c,(r,/ 

nia [ \ nib) 

Cpif,f) = -47r/"/ + r;H'; + ffH', 

Ctifj") = -stt/"/ + /; [2Gf,, + 2g;;, + -G'„ - 



-G'I, + 2Gf,_ + 2G^, 



(28) 



(29) 



6. Numerical Examples 

In order to test the convergence of the algorithm, it is desirable to compare the results to a 
nontrivial exact solution. For a right-hand side consisting of a radially symmetric Gaussian: 

(47n') ' 

we can compute the exact solution to both the Poisson and the biharmonic equations 

Am = / A^v = / 



12 



as well as to the components of the collision operator Cp, Cb- 



u(p) = 
v(p) = 

Cn = 



R 

47Tp 



R = erf 



2 



R- 



47r3/2 



ER 



8^2^,3 l6^3/2^,5/2p' 



Ch 



ER 



Inh'^ 47r3/2^,5/2p- 



(30) 



After a change of variables to cylindrical coordinates, we can find explicit analytic formulas for 
all quantities produced by our solvers (although some of the formulas need to be treated carefully 
to avoid catastrophic cancellations in their numerical evaluation). 
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Example 1: Let us first consider tlie convergence of the solver for a single Gaussian of variance 
V = 0.223. Using the Chebyshev (spectral integration) solvers, note that after an initially rapid 
convergence, the higher derivatives start to diverge due to the ill-conditioning of the linear system 
(FigD. 



Poisson solution convergence - Chebyshev solver - fixed Gaussian variance 




N 

r 



Biharmonic solution convergence - Chebyshev solver - fixed Gaussian variance 




Figure 1: The relative Loo errors for tlie spectral integration based solvers, when applied to the Poisson and biharmonic 
equations (top and bottom, respectively). 
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Using the Green's function based ODE solver, we see that the the solution and its derivatives 
converge spectrally without the loss of precision for fine grids that affects the spectral integration 
based schemes (Fig|2]i. 



Poisson solution convergence ~ Green's function soiver - fixed Gaussian variance 





Figure 2: The relative L„ errors for the Green's function based solvers, when applied to the Poisson and biharmonic 
equations (top and bottom, respectively). 
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Fig. [3]depicts the runtime performance of the solver In the present implementation, the spectral 
integration based code for the Poisson equation requires about 25 seconds for 5 million grid 
points, while the biharmonic solver is about three times slower. This is within a small factor of 
the performance of FFT-based codes for doubly or triply periodic constant coefficient PDEs on 
regular grids. The Green's function based solvers are a bit slower at present. We estimate that 
straightforward optimization/precomputation could result in a factor of 2-3 speed-up. 



Overview of runtimes of all solvers - Good FFT size data points 

300 I 1 1 1 1 




2 4 6 8 10 12 

Number of grid points ^ ^ qG 



Figure 3: Runtime for the Green's function and spectral integration based solvers, when applied to the Poisson and 
bihamionic equations. 
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Example 2: The Coulomb collision operator appears as a source term in the Boltzmann-Fokker- 
Planck equation ([T]i. Not considering convective and electromagnetic terms, it affects the evolu- 
tion of ionized gas consisting of a single species as 

d,f + (convective and electromagnetic terms) - C{f, f) 

For a single species of ion, the Maxwellian distribution (a Gaussian centered at the origin) 
is an equilibrium state. From the formulas (|28^ and d30t . it can be verified analytically that 
C(f'^,f'^) = 0. Computing this result numerically on a 192 x 128 grid we obtain zero to about 
14 digits. Note that Cp and Ci, do not vanish independently. There is a real cancellation between 
the two contributions when inserted into the formula (|28). 



\\Cp(f,f)\\oo 




l|C(/^,/'')l|oc 


3.403 ■ 10-3 


1.361 ■ 10-2 


5.851 ■ 10-'4 



To illustrate the diffusive nature of the collision operator, let us consider a perturbation to the 
equilibrium solution, by constructing an anisotropic Gaussian density, with slightly different 
variances in the r and z directions: 



where Vr - 1.107 and - 1.353. 



1 



Anisotropic Maxwellian 



Difference between the ar 



Anisotropic Maxweliian - Coilision operator 





Figure 4: Anisotropic Maxwellian, its difference from a Maxwellian and its computed collision operator 



Notice that where the anisotropic distribution is too small (the valley in the central plot), the col- 
lision operator is positive, thus it tends to increase /. Where the anisotropic distribution function 
is too large (the peaks in the central plot), the collision operator is negative, tending to decrease /. 
Thus, the collision operator indeed has the effect of moving an anisotropic Maxwellian towards 
an isotropic equilibrium distribution. 
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To demonstrate the full three dimensional solver, we placed three Gaussian source densities at 
different locations in R^, two with positive weight and one with negative weight. The following 
plot shows the solution overlaid on the r - 6 - z grid as a three dimensional contour surface plot 
with an octant "cut out": 




We used Gaussian variances of 0.2, 0.6 and 0.3 centered at the points with Cartesian coordinates 
(4.3,1.2,3 .6), (- 1 . 1 , 4. 1 , -0.8), and (5 .3, 3 .5, -3.2) and weights - 1 .0, - 1 .3, and 1 .2, respectively. 
With a grid in (r, 0,z) of 192 x 64 x 96 on [0, 16] x [0, 27r] x [-16, 16], we obtained 12 digits of 
accuracy in the solution and its gradient. The total execution time was 18 sees, on a single core 
of a 2.5GHZ CPU, using an oversampling factor in z of 4. 
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7. Conclusion 



This paper describes a new fast solver for separable elliptic partial differential equations in cylin- 
drical coordinates that is both fast and high-order accurate, with solution times comparable to a 
few FFTs using the same number of degrees of freedom. Combined with the Rosenbluth for- 
malism, it permits the rapid evaluation of the Coulomb collision operator in kinetic models of 
ionized gases. 

Our solver is particularly useful when the number of azimuthal modes is small. For full three- 
dimensional problems, it is quite likely that Cartesian-based methods will be more effective, 
particularly since one can use fast multipole-based, fully adaptive solvers. Here, we require 
regular grids in the 9 and z directions, which is sufficient for many applications. 
Several open problems remain. One involves the construction of fast, fully implicit collision 
operators, so that large time steps can be taken in the Boltzmann-Fokker-Planck equation (for 
which there is already a significant literature). Another involves the development of methods for 
solving elliptic partial differential equations in complicated axisymmetric geometries rather than 
in free space - that is, interior or exterior to a surface of revolution. These problems are currently 
being investigated, with progress to be reported at a later date. 
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